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MACHINE  SOLUTIONS  OF  PARTIAL  DIFFERENTIAL 
EQUATIONS  IN  THE  NUMERICALLY  GENERATED 
COORDINATE  SYSTEMS 

by 

Z.  U.  A.  Warsi  and  J.  F.  Thompson,  Jr. 

Abstract 

This  paper  presents  a study  in  depth  on  the  numerically  generated  body- 
fitted  coordinate  systems.  All  fundamental  ideas  have  been  developed  from 
a very  basic  point  of  view  which  leads  to  a well  structured  method  of  coordi- 
nate system  control  in  any  region  of  interest.  It  is  also  shown  how  the 
developed  concepts  are  instrumental  in  obtaining  the  machine  solutions  of 
the  field  differential  equations  in  general  dependent  variables  on  the 
generated  coordinate  systems. 

1.  Introduction 

The  accuracy  of  machine  solutions  of  the  partial  differential  equations 
of  mathematical  physics  depends  very  strongly  on  our  ability  to  impose  the 
required  boundary  conditions  as  accurately  as  possible.  It  is  quite  natural 
to  envisage  a coordinate  system  in  which  one  coordinate  curve  passes  exactly 
through  the  points  at  which  the  boundary  values  of  the  field  equations  have 
been  prescribed  in  the  primary  Cartesian  system  (x,y).  Let  the  boundary 
data  curve  be  denoted  as  n=constant.  To  complete  a two-dimensional  coordinate 
system  we  can  orient  the  other  coordinate  in  any  desired  fashion  and  call  it 
the  curve  ^“constant.  With  this  done,  the  first  step  is  to  generate  the 
coordinate  curves  so  as  to  cover  the  entire  field  in  which  the  solution  of 
the  partial  dlfferentail  equation  is  sought.  The  resulting  complexity  of  the 
field  equations  in  these  coordinates  is  relatively  unimportant. 

In  fluid  mechanics,  and  also  in  many  other  areas  of  mechanics  and  physics, 


the  boundary  values  are  usually  prescribed  on  closed  curves.  For  example, 


the  no-slip  condition  of  viscous  flow  is  prescribed  on  a closed  curve  which 
is  the  body  contour  itself.  Thus  the  coordinate  n=constant  is  chosen  as 
the  body  contour  on  which  the  variations  of  x and  y are  known  in  advance. 
Without  loss  of  generality  we  can  enclose  the  body  by  another  closed  curve 
on  which  the  field  values  are  either  prescribed  or  can  be  calculated  quite 
easily.  For  example  in  fluid  mechanics  we  can  enclose  the  body,  past  which 
the  flow  takes  place,  by  another  curve  sufficiently  far  away  from  the  body, 
where  the  free  stream  values  can  be  prescirbed.  The  idea  of  numerical  coor- 
dinate generation  is  to  fill  the  region  between  the  body  and  the  external 
boundary  curves  by  coordinate  lines  in  the  physical  (x,y)  space.  Because 
of  the  closed  region,  it  is  convenient  to  solve  a system  of  elliptic  equations 
of  the  simplest  type,  which  are  the  Laplace  equations. 

The  preceding  ideas  in  one  form  or  another  have  been  used  by  Winslow  [1], 
Barfield  [2],  Chu  [3],  Amsden  and  Hirt  [4],  and  Godunov  and  Prokopov  [5]. 
However,  the  whole  concept  has  been  used  in  a much  organized  manner  to  provide 
a series  of  solutions  in  fluid  mechanics  by  Thompson,  Thames,  and  Mastin  [6], 
[7],  [8],  [9],  [10],  [11].  References  [7]  and  [11]  contain  a thorough  discus- 
sion on  the  actual  computational  aspects  and  computer  programming  methods  both 
for  the  coordinate  generation  and  the  numerical  solutions  of  the  Navier-Stokes 
equations. 

In  this  paper  we  summarize  the  fundamental  ideas  of  the  method  of  coor- 
dinate generation  which  naturally  lead  to  a well  structured  method  of  the 
coordinate  system  control  in  the  desired  regions  of  the  field.  Further  the 
use  of  numerically  generated  coordinates  in  the  machine  solutions  of  the 
Navier-Stokes  equations  written  in  the  contravariant  components  of  the  velocity 
vector  is  demonstrated. 

2.  Method  of  Numerical  Coordinate  Generation 

Let  C ” £(x,y)  and  n ”n(x,y)  be  two  continuously  differentiable  functions 


i 


2 


of  the  Cartesian  coordinates  (x,y).  The  functions 
£(x,y)  = , n(x,y)  = C2 

define  curvilinear  coordinates  in  the  xy-plane.  By  varying  the  constants 
and  C2  we  can  cover  the  entire  xy-plane  by  the  curvilinear  nets  or  meshes. 

The  choice  of  the  functions  £ and  n is  quite  arbitrary  because  the  final  out- 
come of  the  solution  of  any  field  differential  equation  either  in  the  (x,y) 
or  (£,n)  system  remains  the  same.  However,  the  choice  of  a suitable  coordinate 
system  is  usually  governed  by  the  geometrical  configuration  of  the  physical 
field  particularly  with  a view  on  satisfying  the  boundary  conditions  in  as 
simple  and  exact  manner  as  possible.  Since  most  of  our  discussion  in  this 
paper  is  in  the  context  of  the  viscous  flows  past  finite  bodies,  the  most 
important  boundary  contition  to  be  satisfied  by  the  Navier-Stokes  equations 
is  the  vanishing  fluid  velocity  at  a stationary  finite  body. 

For  the  convenience  of  satisfying  the  exact  boundary  conditions  on  the 

body  surface,  we  choose  one  coordinate  line  n=n  =constant  as  the  body  contour. 

o 

If  the  body  is  finite  then  n=n  is  a closed  curve.  We  now  enclose  the  body 

o 

by  another  closed  curve  and  call  it  the  outer  boundary,  where  numerically 

n >q  . The  region  q < q<q  must  now  be  filled  by  coordinate  curves  by  some 

CX>  Q 0 — — OO 

method.  Because  of  the  closed  region  under  consideration  it  is  natural  to 
specify  the  determining  differential  equations  for  £ and  n as  elliptic  equa- 
tions to  be  solved  under  proper  boundary  conditions  at  the  body  and  at  the 
external  boundary.  The  simplest  elliptic  equation  is  the  Laplace  equation. 

With  this  in  mind  we  pose  the  problem  of  solving  the  Laplace  equations  for 


£ and  n with  x and  y as  independent  variables  under  the  Dlrichlet  boundary 
conditions.  Let  be  the  curve  defining  the  body  contour  q=qQ  and  T2  be  the 
curve  defining  the  outer  boundary  q«q  in  the  xy-plane.  Fig.  1.  The  elliptic 





boundary  value  problem  is  then 

v2e  = o , (i) 

v2n  = 0 , (2) 

£ = fQ(x,y)  , n = n0  on  1^  , (3) 

5 = ^(x.y)  , n = on  r2  , (4) 

To  fix  ideas  we  consider  the  example  of  plane  polar  coordinates  (r,9), 
where  it  is  known  that 

x = r cose  , y = r sin0.  (5) 

It  is  easy  to  show  that  V20  = 0 but  V2r^  0.  However,  the  variables  defined 


e - e , 


n = Jtnr  , 


where 


x = r cos 0 = e*'nr  cose  = encos£;  , 


y = r sin0  = e^nr  sine  = ensin£ 


are  such  that  V2£  = 0 and  V2n  = 0.  If  we  take  the  inner  and  outer  bounding 

circles  of  radii  r^  and  r 2 respectively,  then 

nQ  = = «.nr2.  (9) 

In  this  case  we  also  have 

f = f tan  ^ ^ . 
o » x 

Thus  the  solution  of  the  elliptic  problem  (l)-(4)  for  the  circular  annulus  is 
£ = tan-1  ^ , n = -|«,n(x2  + y2) 


The  above  example  illustrates  that  it  should  be  quite  convenient  to 
solve  Eqs.  (l)-(4)  numerically  in  those  cases  in  which  it  is  possible  to 
specify  the  values  of  nQ  and  by  simple  analytic  means.  Unfortunately  the 
list  of  such  body  and  outer  boundary  shapes  and  configurations  is  painfully 


short.  Thus  the  method  of  numerical  coordinate  generation  for  arbitrary 
shaped  bodies  and  outer  boundaries  through  Eqs.  ( 1) — (4 ) will  not  be  very  prac- 
tical. The  source  of  difficulty  lies  not  so  much  in  the  specifications  of  fQ 
and  f^  but  in  the  specifications  of  nQ  and 

To  overcome  these  difficulties,  and  also  for  other  important  purposes  to 
be  described  later,  we  interchange  the  roles  of  the  dependent  and  independent 
variables  in  Eqs.  (l)-(4).  Thus  using  the  relations 

= Jy  xn  =-jy  yc = _Jv  \ = Jy  (10) 

where 

J - /i  - x6yn  - *nyt  , (u) 

and  the  implicit  differentiation  formulae 

'k=^'k+\J^  etC*  ’ (12) 

we  obtain 

g22  ~ 2g12  x£n  + ell  xnn  0 ’ {13) 

822  ' 2g12  + 811  ynn  ‘ ° ' (U) 

(Equations  (13)  and  (14)  can  also  be  written  down  directly  by  first  writing 

and  then  4>=  n in  Eq.  (Al-33)  and  then  solving  them  simultaneously).  For  the 
definitions  of  the  quantities  g^  refer  to  Appendix  A-l. 

The  transformation  (10)- (12)  implies  that  the  body  and  the  external 
boundary  contours  in  the  xy-plane  have  been  mapped  on  to  the  £n-plane.  In 
other  words  we  can  say  that  the  contours  in  the  xy-plane  have  been  opened 

up  to  form  the  straight  lines  n=n  =const.,  and  n=n  =const.  in  the  Cn-plane. 

0 00 

This  can  be  achieved  by  imagining  a cut  connecting  the  body  and  the  external 
boundary  in  the  xy-plane  as  shown  in  Fig.  1,  such  that  all  functions  and  their 
derivatives  are  continuous  in  crossing  the  cut. 


The  boundary  conditions  for  Eqs.  (13)  and  (14)  are 


5 


(15) 


x = F..(£,n  ),  y = F0(£,n  ) on  T * , 


X = G1(C,noo),  y = G2U,r\m)  on  , 


(16) 


where  as  shown  in  Fig.  2,  T • and  are  the  images  of  the  body  and  the 


external  boundary  contours  in  the  £n-plane.  Obviously  no  boundary  condi- 

* , r,  V 


tions  can  be  specified  on  and  . 

Equations  (13)  and  (14)  are  quasilinear  uniformly  elliptic  partial 
differential  equations,  as  proved  in  Eq.  (A2-6)  of  Appendix-A2. 


The  appearance  of  rio  and  in  (15  and  (16)  is  now  purely  symbolic. 


denoting  the  names  of  the  body  and  the  external  boundary  respectively.  Given 
the  body  and  the  external  boundary  contours,  we  can  always  establish  the  val- 
ues of  x and  y either  graphically  or  analytically  for  any  desired  distribution 
of  the  £-values.  The  n-values  can  also  be  chosen  arbitrarily  to  form  rectangu- 
lar meshes  in  the  £n-plane. 

A question  which  naturally  arises  at  this  stage  is  this:  how  to  choose 

the  variations  of  £ and  n or  how  to  label  them  in  covering  the  whole  £n-plane? 
An  answer  to  this  question  is  that  the  whole  operation  can  be  more  practically 
oriented  if  we  can  establish  a correspondence  between  the  actual  curvilinear 
coordinate  values  (£,n)  and  the  set  of  integers.  To  effect  this  transition, 
we  introduce  new  variables  (x,o)  through  a linear  transformation 


£ = ay  + £q  , n = bo  + , 


(17) 


where  a and  b are  constants.  Carrying  out  the  transformation  (17)  in  (13) 
and  (14),  we  have 


822  Xxx  " 2g12  XXo  + 811  Xoo  = ° ’ 


822  yXX  ‘ 2i12  yX0  + 811  yoo  = 0 * 


(18) 

(19) 


where  g^  are  functions  of  (x.o),  viz., 


2 - 


811  * Xx  + yX  ’g12  = Vo  + Vo’  822  = xo  + V 


For  these  and  other  definitions  refer  to  Appendix-Al.  Here  it  must  be 
noted  that  for  any  function  F(£,n), 


& 


F = a(F  ) 

X e C=ax+C, 


F = b(F  ) 

o n n=bo+n 


(20a) 


FXX  " 3 (FCC)C=aX+C0’ 


Foo  ^ ^Fnn^n=bo+n  ’ 

O 


F = ab  f F ], 

xa  £n  C=ax+C, 


(20b) 


n=bo+n 


Equations  (18)  and  (19)  are  exactly  of  the  same  form  as  Eqs.  (13)  and 
(14).  However  a great  simplification  can  be  achieved  if  x and  o are  treated 
as  integers,  for  then  irrespective  of  the  numerical  values  of  a and  b we  can 
proceed  either  in  the  £ or  n directions  through  consecutive  integers. 

Further,  according  to  the  right  hand  sides  of  Eqs.  (20  a,b)  the  first  deriva- 
tives appearing  in  Eqs.  (18)  and  (19)  can  be  approximated  by  forward,  back- 
ward, or  central  differences  without  specifying  a or  b.  Similarly,  one  can 
approximate  the  second  derivatives  in  Eqs.  (18)  and  (19)  by  a central  dif- 
ference. All  higher  derivatives  can  similarly  be  approximated  by  the  proper 
differences.  Under  no  circumstance  should  the  quantities  like  etc.  be 
interpreted  as  derivatives  with  respect  to  integers! 


The  machine  solution  of  Eqs.  (18)  and  (19)  can  now  be  easily  performed 
by  employing  the  proper  differences  for  the  derivatives  without  specifying 
the  step  sizes.  The  same  argument  is  of  course  true  for  Eqs.  (13)  and  (14) 
because  they  are  also  homogeneous.  However,  it  must  be  emphasized  that  the 
approach  taken  here  is  fundamental  in  building  the  essential  concepts  and 
establishing  the  versatality  of  the  method.  The  reason  for  developing  the 


1 


1 


j 


I 

I 


1 
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concepts  this  way,  rather  than  the  one  originally  proposed  , is  that  the  method 
of  numerical  coordinate  generation  and  the  solutions  of  other  differential 
equations  (field  equations),  either  homogeneous  or  inhomogeneous  over  the 
generated  meshes,  must  be  considered  together  and  not  separately.  It  will 
be  shown  later  that  the  same  concepts  are  applicable  when  the  step  sizes  in 
the  Cn-plane  are  variable  and  not  fixed  like  a and  b. 

The  main  utility  of  numerically  generated  body-fitted  coordinates  actually 
lies  in  the  availability  of  meshes  cr  nets  in  the  £r|-plane  on  which  any  dif- 
ferential or  integral  equation  representing  a field  can  be  solved.  First 
the  differential  or  integral  equation  must  be  transformed  to  a general 
coordinate  system  (£,n)  and  then  tc  (v.a).  If  the  equation  in  the  (y.o) 
variables  does  not  show  any  explicit  dependence  on  a and  b then  we  can  make 
a finite  difference  approximation  of  the  equations  without  specifying  the 
step  sizes.  This  statement  is  equally  true  for  the  variable  mesh  sizes. 

Later  in  the  paper  we  have  demonstrated  that  in  the  case  of  the  Navier-Stokes 
equations  the  dependent  variables  can  always  be  selected  in  such  a way  that 
the  step  sizes  do  not  appear  explicitly.  The  most  important  point  to  mention 
is  that  in  the  midst  of  these  transformations  and  adjustments  of  terms,  the 
physical  outcome  of  the  solution  is  not  disturbed. 

t 

In  their  original  proposal  the  authors  [6],  [7],  treated  the  actual 
curvilinear  coordinates  (£,n)  as  integers.  The  main  reason  for  not  intro- 
ducing the  above  concept  was  that  the  Navier-Stokes  equations  which  they 
solved  had  the  Cartesian  components  of  velocity  as  independent  variables 
so  that  the  step  sizes  a and  b did  not  appear  in  their  equations. 


1 


a 


■ 
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2.1  Variable  Meshes  and  Coordinate  Contraction: 


i 


To  see  the  effect  of  the  step  size  more  clearly  and  to  have  a method 
of  changing  the  step  sizes,  we  make  a general  transformation  in  the  basic 
elliptic  system  Eqs.  (13)  and  (14)  of  the  form 


C = f-i  (x)  + £ 

i o 


f0(o)  + T) 


where 

and 


C = Co  at  X 


Writing 


o 

df. 


0 , 
at  o = 0 . 


df. 


A(x)  = 


dX 


e(o)  = 


do 


^ _ 1 

d£  A 


do 

dn 


so  that  the  derivatives  transform  as 


x = x /A  , 

£ X 


x = x / 0 , 

n a 


J = J/9A  ; J = x y - x y 
X o a7x 


Further  noting  that  A = A(y)  and  0 = 0(o),  we  have 

Ax  _ 

- = (x  _ X -X)  m 2 

1 XX  A )/X  * 


cc 


xr  = x /0A  , 

tn  e x , 

x = (x %-^)/02 

nn  oo  0 


(21) 

(22) 


(23a) 


(23b) 


and  similar  expressions  for  the  second  derivatives  of  y. 


ye 

r< 

X 

II 

yn 

= ya/e  , 

11 

- hin2 

; ®11 

= x2  + 

X 

2 

yx  * 

(24) 

12 

• *12/61 

’ ®12 

= X X 

X o 

+ yxyo  • 

’22 

gz2/0 

g22 

= x2  + 
a 

2 

yo  * 

g 

= I/02a2 

: i = 

iiii22 

- <i12)2  * 

Substituting  (24)  and  (25)  in  (13)  and  (14),  we  have 


*22xxx  - 2i12  XXo  + 5U  x«,o  ‘ Pxx  + QX»  ’ 


i22yxx  ' 2*12  V + *11  y™  ' Pyx  + ^0  ’ 


where 


®22  811 
p = irAx  ’ Q = 


(26) 

(27) 

(28) 


It  is  a matter  of  direct  verification  that  the  differential  equations 
V2X  =-R(x,o)  , (29) 


and 


V2o  =-S(x,o)  , (30) 

where 

V2  = Jl  + R * (x2  + xz)A  A , S = (o2  + o2)0n/e 

lx7  dy7  , X Ay  X x y a 

transform  to  the  same  forms  as  Eqs.  (26)  and  (27).  Differential  equations 

of  the  form  (29)  and  (30)  yield  solutions  which  are  called  superharmonic  if 
R>0  and  S>0.  As  is  proved  in  Appendix-A3,  the  superharmonic  functions  like 
harmonic  functions  always  satisfy  a minimum  principle  and  are  bounded  from 
below  by  the  values  of  the  harmonic  functions  satisfying  the  same  boundary 
conditions . 

As  in  the  case  of  uniform  step  sizes,  the  finite  difference  approximation 
of  Eqs.  (26)  or  (27)  also  does  not  depend  on  the  step  sizes.  Further,  because 
of  the  vai tables  X and  0 we  can  now  control  the  mesh  spacings  both  in  the 
£ and  n directions  merely  by  specifying  the  functions  f^(x)  and  f2(o).  The 
functions  f^  and  f^  must  be  specified  as  continuous  and  preferably  analytic 
functions  of  their  arguments.*  Care  should  be  taken  in  not  interpreting 


t 

— oo<  ^ <oo 


It  is  preferable  to  have  both  f^(x)  and  f2(°)  analytic  in  the  range 
and  -<=<a<oo  to  admit  all  integral  values  for  x and  o. 
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X and  a as  integers  till  all  their  needed  derivatives  have  been  obtained 
in  the  usual  way.  Once  A^  and  0^  have  been  obtained  and  inserted  in  equations 
(28),  the  variables  x and  a can  then  be  interpreted  as  integers  in  the  solu- 
tions of  Eqs.  (26)  and  (27).  The  following  example  illustrates  this  point. 

Suppose  we  need  to  compact  the  mesh  in  the  neighborhood  of  the  body 
surface  to  have  a better  resolution  of  the  viscous  effects  in  a fluid  flow 
problem.  Since  no  contraction  is  needed  along  the  ^-direction,  we 


f1(x)  = ax  , = 0, 


where  a = constant. 


For  the  function  f^(o)  let  us  prescribe  an  analytic  function 
f2(o)  = bo<°  * , 

where  k and  b are  constants.  With  (31)  and  (32)  we  have 

A = a = const.  , 

9(o)  = b<°  ^(l+allnK)  , 

0 = bK°  1 (2+o£n<) UnK . 

o 


Therefore  from  (28) , 


P = 0 , 


g11[2«-nic  + a(«.n<)  ] . 

Q = 1+afcnK 

It  is  interesting  to  note  that  the  constant  b does  not  appear  in  Q.  Insert- 
ing (35)  in  Eqs.  (26)  and  (27)  we  get  the  equations  in  which,  by  approximating 
the  derivatives  by  differences,  we  can  treat  both  x and  o as  integers.  Equa- 
tions (26)  and  (27)  along  with  the  values  of  P and  Q given  in  (35)  have  been 
solved  for  the  case  of  an  ellipse  with  tc =1 . 1 and  the  generated  coordinates 

4. 

have  been  shown  in  Fig.  3.  The  generated  values  of  x and  y for  this  case  match' 
quite  closely  with  the  available  exact  conformal  solution. 


Note  that  we  could  match  the  generated  solution  with  available  exact 
solution  because  given  the  semi-major  and  minor  axes  of  both  contours  we  could 
find  the  constant  b. 
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3.  Application  to  the  Navier-Stokes  Equations 

In  this  section  we  demonstrate  that  in  the  case  of  the  Navier-Stokes 
equations  the  form  of  the  equations  whether  written  in  the  actual  curvilinear 
coordinates  or  in  the  generated  curvilinear  coordinates  remains  the  same. 

This  result  implies  that  in  the  solution  of  the  Navier-Stokes  equations  the 
step  sizes,  whether  constant  or  variable,  are  not  to  be  specified  from  out- 
side. The  effect  of  the  step  size  is  however  built  into  the  generated  coor- 
dinate meshes  which  have  already  been  obtained  a' priori  by  solving  Eqs . (26) 
and  (27). 

The  Navier-Stokes  equations  in  non-dimensional  form  written  in  general 
curvilinear  coordinates  and  with  the  contravariant  components  of  the  velo- 
city vector  (v1)  and  1/e  as  the  Reynolds  number  are 

|r  + (>/gpvi)  = 0 , (36) 

9 , i,  .1  3 / /-  i js^Pijk 

(pv  ) + j (p/gv  vJ)  + T pvJv 

9xJ  3 

1 9 , j-  ij,  „i  jk 

= - --  — r (p*^  8 ) - r Pg 

✓g  Jk 

+ e [— (/gsg1^)  + — “t  (^  di3) 

4 K5  3Cj 

+ sr^  gjk  + djk)  , (37) 

where  p is  the  density  (absolute  scalar), 

s = p'div  v , 

d1J  - ,<8lkv^  + v‘l  ) . 

u and  u'  being  the  first  and  second  coefficients  of  viscosity,  respectively. 
For  all  other  notation  and  definitions  refer  to  Appendix-Al. 

In  two  dimensions,  Eq . (37)  represents  two  equations  (i»l,2).  Intro- 
ducing the  transformation  (21)  and  using  (24)  it  can  be  shown  that 
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r 


L I 

r 


r Ur  1-^-x 

11  A 11  x2  X ’ 


r 2=  i f 2_  J_ 

22  6 22  .2  o , 

U 


r 1-±r  1 

22  0 22 


r 2 = — f 2 

11  A2  11 


r 1. 

If  1 

ri2 

0 12 

9 

2 

ri2  = 

I f 2 

A ri2 

> 

r.  1 . 

r.  2 

1 - 

F11 

r = 

12 

2Ag  8X 

„ 1 . 

_ 2 

1 - 

ri2  + 

r22  = 

20g  8o 

(38) 


0 


where  T 1 has  the  same  form  In  (x,o)  as  r*  has  in  (£,n). 


If  we  now  define  new  contravariant  components 


1 v1  2 v2 

u = T * u = T 


(39) 


then  by  substituting  (38)  and  (39)  in  (36)  and  (37)  the  equations  retain 
the  original  form  except  that  the  independent  variables  are  (x,o)  and  the  g's 
and  r's  are  over  barred.  At  no  place  A and  0 appear  explicitly  in  the  equa- 
tions. Therefore  if  we  treat  x and  a as  integers  then  all  the  derivatives 
can  be  approximated  by  their  corresponding  differences. 

It  was  mentioned  earlier  that  by  using  all  these  transformations  the  phys- 
ical outcome  of  the  solution  is  not  disturbed.  The  following  example  is 
sufficient  to  demonstrate  this  fact. 


1 2 

If  v and  v are  the  contravariant  components  of  the  velocity  vector  v 


then  the  magnitude  of  the  velocity  is  given  by 


MZ  “ 8i  i (v1)2  + 2g  vV  + g00(vU  . 


1 2 


2.2 


’ll'*  ' ' '"612  " v ' b22' 

Using  (2A)  and  (39),  the  same  magnitude  is  now  given  by 


(40) 


lv|2  ” 811(ul)  + 2g12  u1u2  + g22(u2)2 


(41) 
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| 3 


It  is  a matter  of  direct  verification  that  the  Navier-Stokes  equations 
written  in  any  other  form,  such  as  in  the  covariant  or  cartesian  components 
of  the  velocity  vector,  can  be  solved  in  the  generated  coordinate  system 
without  any  specification  of  the  step  sizes.  In  passing  we  may  note  that  in 
two  dimensions  the  relations  between  the  Cartesian  and  the  contravariant 
components  of  the  velocity  vector  v are  simply 

„ 1^2  1^2 

U-XV  + X V = X u + X u , 

Z n X o 

„ 1 t 2 1 . 2 

V = y^v  + y^v  = y^u  + y0u  . (42) 

4.  Truncation  Errors 

In  any  finite  difference  approximation  it  is  imperative  to  have  an  estimate 
of  the  truncation  errors  to  establish  the  accuracy  and  consistency  of  the 
difference  scheme.  To  investigate  the  problem  of  truncation  errors  in  the 
present  case  we  consider  a function  F (£,n)  of  the  coordinates  (£,n).  For 
simplicity  consider  that  £ is  fixed  and  n assumes  consecutive  values  r^,  r^. 

Let  us  use  the  following  notation 


F.  - FU,  n.)  , Orf) 
i i 3n  n=n. 


, etc.,  i = 1,2,3, 


Using  Taylor’s  expansion,  we  have 

F1  ■ F2  - <VW  + 2'VV2  F2"  - I'VVV"  + • 

F3  - F2  + fyV'i'^YV2  f2"  + t<vn2,3Y" + • 


-1— i - (il)  . i(„  _2„  +n  )F  --  [(n3-n2)  +(n2-n,)3)] 

n,-n,  2Vnl  ZVn3;t2  + =— ± F2  + 

3 1 6(n,-n,) 


Equation  (43)  shows  the  difference  between  the  computed  and  true  derivative, 
viz.,  the  truncation  or  the  discretization  error  in  the  first  derivative,  in 
terms  of  the  consecutive  values  n^»  n ^ and  n^.  To  find  the  order  of  magnitude 
of  the  right  hand  side  in  Eq.  (43),  we  introduce  the  transformation 
n - f-(o)  + n , 

L O 
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Using  Taylor's  expansion  of  (44)  about  o , we  get 


and  a subscript  2 on  parentheses  denotes  the  value  at  o=o 


must  be  emphasized  again  that  6(a)  and  all  its  higher  derivatives  are  obtained 


by  differentiation  of  the  continuous  function  f (a) 


If  we  now  interpret  o as  assuming  consecutive  integral  values  then 


To  have  a numerical  estimate  of  the  terms  on  the  right  hand  sides  of  (45) 


and  (46)  we  have  to  prescribe  an  f_(o)  to  determine  9 and  its  derivatives.  For 


example  the  function  f,(o)  as  prescribed  in  (32)  shows  that  both  9 and  9 are 


of  the  order  of  b,  which  is  implicitly  assumed  to  be  a very  small  number 


The  above  analysis  also  points  out  that  the  choice  of  the  functions  f 


and  f.  must  be  made  with  the  following  two  criteria 


1.  Both  f.  and  f must  be  analytic 


The  values  of  9 and  all  its  derivatives  be  small  for  all  values  of 
X and  o. 


5.  Non-Steady  Coordinates 


Body-fitted  coordinates  are  subject  to  change  with  time  if  the  body 
and/or  the  external  boundary  contours  change  with  time.  The  physical  situa 


tions  in  which  the  non-steady  nature  of  the  coordinates  must  be  considered 
are,  for  example,  when  coordinates  are  attached  to  a pulsating  body,  to  a 
propagating  blast  or  shock  wave,  free  surface,  etc. 

It  will  be  shown  in  this  section  that  by  using  implicit  differentiation 
rules  the  non-steady  nature  of  the  coordinate  system  can  be  brought  into  the 
time  derivative  terms  of  the  differential  equations  which  are  to  be  solved  on 
the  generated  coordinate  system.  For  each  time  step  the  coordinates  have  to 
be  generated  or  upgraded  by  using  the  same  equations,  viz.  (26)  and  (27)  under 
the  changed  boundary  data.  The  greatest  facility  with  this  approach  is  that 
no  interpolations  are  needed  between  two  successive  time  intervals  and  the 
numerical  solution  proceeds  with  the  same  algorithm  as  when  the  coordinates 
are  not  changing  with  time. 

Referring  to  Fig.  1,  let  the  point  p away  from  the  body  has  Cartesian 
coordinates  (x,y).  In  fluid  flow  problems  (x,y)  are  the  Eulerian  variables 
and  represent  a fixed  point  in  the  physical  space.  Thus,  no  matter  how  the 
body  moves,  or  the  oncoming  flow  takes  place,  the  total  time  derivative  of 
x and  y are  zero.  That  is 


3T-°  • ■ <47> 

Let  f(x,y,t)  be  a flow  quantity.  Then  in  the  body-fitted  coordinates 
if  the  body  and  the  external  boundary  are  kept  fixed,  then  no  matter  how  the 
flow  takes  place,  the  partial  derivative  of  f is  given  by 

(— ) « (M) 

V3t'x,y  v9t'ii,n  (48) 

since  in  principle  x » x(£,n),  and  y « y(£,n).  On  the  other  hand  if  the  body- 

fitted  coordinates  are  time  dependent,  then 

K “£  (x,y,t),  n = n(x,y , t)  , (49) 

from  which  we  deduce  that 

(|f)  „ - (If)  + M.  + fn,  , (50) 

at  x,y  at  £,n  S t n t 

where  as  before  the  variable  subscripts  denote  the  partial  derivatives. 


Now  from  (49),  we  also  have  the  inversion 
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x * x(£,n,t)  , y = y(£,n,t)  . 

Thus  using  (47),  we  have 


(51) 


xt  + Vt  + xnnt  = 0 * (52a) 

yt  + Vt  + Vt  = 0 ’ (52b) 

which  on  solution  give 

£t  * (ytxn  - x^^/J  » (53a) 

nt  = (xtyc  " ytxC)/J  ' (53b) 

Substituting  (53)  in  (50),  we  have 

(JiK,y  " #C,n  + J<V{  ' VnK 

+ J<fs%  ' fnVyt  ’ (54) 

Equation  (54)  can  now  be  used  to  replace  the  partial  time  derivatives  in 
the  Navier-Stokes  equations  written  in  the  Cartesian  velocity  components  (U,V), 
or  in  the  stream  function  i|».  The  same  equation  can  also  be  used  to  replace 
the  time  derivatives  of  density  and  energy  in  the  compressible  flow  equations. 
However,  some  analysis  is  required  to  replace  the  time  derivatives  appearing 
in  the  Navier-Stokes  equations  when  written  in  the  covariant  or  contravariant 
components  of  the  velocity  vector. 

gv 

Let  v be  the  velocity  vector  so  that  ~ is  the  partial  time  derivative 

O L 

term  in  the  vector  form  of  the  Navier-Stokes  equation.  If  time  dependent 
coordinates  are  used  then  the  form  of  the  time  derivative  becomes 
3v  9v  grk 

(atV  + 7ik  h ’ (55) 

at. 


3V 

and  not  simply  as  (— )^  . Here  for  ease  of  writing  the  expressions  we  have 

-i 

denoted  the  curvilinear  coordinates  by  £J  (j»l,2)  . 


i 

v’k  ®i 
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where  v,*  is  the  covariant  derivative  of  the  contravariant  component  v*,  and 
a1  represent  the  covariant  base  vectors.  From  (Al-16)  and  (Al-18) , we  have 

v*£  - + rki vr  • (57> 

The  vector  v written  in  the  contravariant  components  v is 

v = aiv1  (58) 

Substituting  (58)  in  (55)  we  get  on  using  (56) 

..  i 9 a k 

~dt  ®i  + ~Tt  V + V,k  ~i  • (59> 

0 

Multiplying  each  term  of  (59)  scalarly  by  the  contravariant  base  vector  a 

and  recalling  the  property 

l .1 

3m‘a  = 6 » 

~m  - m 


(59)  becomes 


9V1  . k 3*k 

TT  + V TT 


i . i 9?" 
§ + v.u  TT 


9t  9t  " ’k  9t  ’ 

where  repeated  indices  always  imply  summation. 

1 2 9£^ 

Denoting  £ =£  and  5 =n,  the  values  of  in  (60)  are  those  as  given  in 
Eqs.  (53).  It  must  be  noted  that  if  v is  replaced  by  a scalar  function  then 
the  middle  term  in  (60)  does  not  appear  and  the  expression  coincides  exactly 
with  (54). 


Based  on  the  equations 

i ij 
a =gJaj 


k 3*k  i 

1 3'1  i 2 3'2 

V 9t  ? 

V 9t  5 + v 3t 

for  i = 1: 

k 3§k  1 

1/ 

v 9t 

* v <x? + ’ 

nt  n nt  n 
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For  i = 2: 


k 3*k 


?2  = 


6.  Generalizations  and  Extensions 

For  the  sake  of  clarity  of  exposition  the  preceding  development  has 
been  confined  only  to  doubly  connected  regions,  i.e.,  the  regions  formed  by 
a body  and  an  external  boundary.  The  same  concepts  of  coordinate  generation 
are  applicable  to  mulitply  connected  regions.  Given  any  number  of  bodies, 
we  can  join  the  bodies  through  simple  cuts  and  the  last  body  can  be  joined 
with  the  external  boundary  by  a cut.  The  transformation  from  the  physical 
(x,y)  space  to  the  transformed  (£,n)  space  then  yields  discrete  segments  on 
the  n=const.  line;  each  segment  representing  one  body.  The  computer  algorithm 
for  such  cases  has  been  throughly  developed  by  one  of  us  (JFT)  and  a detailed 
discussion  is  available  in  Ref.  [11]. 

The  most  difficult  and  also  rewarding  area  of  extension  of  the  concepts 
of  coordinate  generation  lies  in  the  three-dimensional  situations.  It  is 
hoped  that  the  ideas  presented  in  this  paper  may  be  of  help  in  clarifying 
some  concepts  in  three  dimensions. 
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Appendix  - A1 


In  this  Appendix  we  summarize  some  results  from  the  theory  of  tensors 
which  are  always  needed  in  the  formulation  of  problems  in  general  curvilinear 
coordinates.  For  further  details  refer  to  Sokolnikoff  [12].  Some  derivations 
can  also  be  found  in  Ref.  [13]. 


Notation  and  Definitions 

Following  the  standard  practice  in  tensor  theory  we  shall  use  variables 
both  in  subscript  and  superscript  index  notation.  In  what  follows,  the 
repeated  indices  always  imply  summation. 

In  the  three-dimensional  Euclidean  space  we  introduce  a Cartesian 
coordinate  system  denoted  as  x.,  where  i=l,2,3.  Embedded  in  this  Cartesian 


space,  we  introduce  a general  curvilinear  coordinate  system  denoted  as  £ , 


where  i=l,2,3.  Thus  x^  are  functions  of  £ and  vice  versa. 


The  magnitude  of  the  displacement  vector  dr  denoted  as  ds  when  referred 


to  the  Cartesian  coordinates  is  given  by 


(ds)2  = 6jdx^dxj 


(Al-1) 


while  referred  to  the  general  curvilinear  coordiantes  it  is  given  by 


(ds)2  = g.jd£1d£J 


(Al-2) 


In  (Al-1)  6^  are  the  Kronecker  deltas  defined  as 


6T  = 1 if  i=j 

j 


= 0 if  ijtj  . 


(Al-3) 


In  (Al-2)  g_  are  the  covariant  components  of  the  metric  tensor.  It  is  a 


symmetric  tensor,  viz.. 


8ij 


(Al-4) 


Equating  (Al-1)  and  (Al-2),  we  can  easily  establish  the  formula 


3xk  3xk 


’ij 


9C1  3£j 


(Al-5) 


ij 


The  contravariant  components  of  the  metric  tensor  are  g J where 


A-l 


(Al-6) 


! 


r i 

.1 


J 


ik  xi 

8 gr  = 6' 


’jk 


Thus 


ij  aci 

8 ~ 3x,  3x 

k k 


(Al-7) 


The  determinant  of  the  matrix  {g  } is  denoted  as  g. 


g = det(gi  ) 


(Al-8) 


The  Jacobian  determinant  of  the  transformation  is 


J = /g  = 


aCxj^,  x2>  x3) 


12  3 

a(r,  c,  o 


(Al-9) 


where  J^O  means  that  the  transformation  is  not  singular. 

In  the  process  of  differentiation  of  vectors  and  tensors  in  general 
coordinate  systems  certain  quantities  known  as  the  Christoffel  symbols 


appear,  which  for  any  values  of  i,j  and  k are  given  by 
i 


r - , + 
? V 


A ’ 


(Al-10) 


3£ 


In  particular 

" k - 1 3 <4)  - -h 


kr 


Si  9Cr 


2g  3£r 


(Al-11) 


The  quantities  are  symmetric  in  the  lower  two  indices. 


r 1 = r 1 

jk  kj 


(Al-12) 


Differentiation  of  Vectors  and  Tensors 
In  general  coordinates  the  subscripts  denote  the  covariant  components 
while  the  superscripts  denote  the  contravariant  components.  In  Cartesian 
coordinates  since  there  is  no  distinction  between  the  covariant  and  contra- 
variant components,  the  indices  have  no  special  meaning.  Thus  the  subscript 
on  x in  all  the  previous  formulae  is  merely  a symbol  and  should  not  be  consid- 
ered a covariant  index.  Below  we  give  some  examples. 

Let  a be  a vector,  then  it  can  either  be  expressed  in  terms  of  its 
covariant  components  a^  or  the  contravariant  components  a*. 

Let  A be  a tensor  of  second  order.  It  can  be  expressed  in  three 

A- 2 


distinct  ways: 


, the  covariant  components, 

A*^  , the  contravariant  components, 

A^  , the  mixed  covariant-contravariant  components  . 

The  displacment  vector  dr  has  components  dx^,  or  (dx^ , dx£,  dx^)  in  the 

Cartesian  coordinates  without  any  significance  attached  to  the  subscripts.  On 

i i 2 3 

the  other  hand  dr  has  the  contravariant  components  d£  , or  (d£  , d£  , d£  ) in 
any  general  coordinate  system  and  actually  this  was  the  reason  for  denoting 


the  general  coordinates  by  a superscript,  though  £ themselves  are  not  the 


contravariant  components.  Thus 

i i 

dr  = — r dr  = a dr 

3£ 


(Al-13) 


The  quantities  a^,  i=l,2,3,  are  called  the  covariant  base  vectors.  Obvi- 


ously 


o — 3 *3 

ij  ~i 


(Al-14) 


Let  <J>  be  a scalar  function,  then 
b = grad  tp 

is  a vector  whose  covariant  components  in  any  coordinate  system  are 

b = J±_  (Al-15) 

1 H1 

Let  v be  a vector,  then  its  gradient  is  a second  order  tensor  denoted  as 
A = grad  v . 


The  mixed  components  of  A are  given  by 

.i  9v*  , _ i r 

*k  ' TT  + fkr  v • 


(Al-16) 


The  covariant  components  of  A are  given  by 
3v 

Alk  ' vr  ' <A1-17) 

The  right  hand  sides  of  (Al-16)  and  (Al-17)  are  called  the  covariant  deriva- 


tives of  the  contravariant  and  covariant  vector  components  respectively. 
The  covariant  derivative  is  always  denoted  by  a coma  (,)  preceding  an 


index.  Thus  in  (Al-16),  (Al-17) 

. i i 
\ = v’k 


Ail  = vi,k 


(Al-18) 


(Al-19) 


Covariant  differentiation  is  also  called  absolute  differentiation  becasue  it 
holds  in  any  coordinate  system  including  the  Cartesian  system  where  the 
Christoffel  symbols  are  identically  zero.  The  divergence  of  a vector  is 
obtained  by  setting  k=i  in  (Al-18)  and  then  invoking  the  summation  convention. 


div  v = v 


(Al-20) 


The  formulae  for  the  covariant  differentiation  of  second  order  tensors  are 


3A., 

a = — — - r r a - r r a 

ik,i  ..i  in  rk  ki  ir  ’ 

a?’  = ^ + r ^ A.r  - r r a^ 

\,i  ,_i  ir  \ ki  r 

+ r,1  Ark  + r k Air  . 

,i  ~ri  ir  ir 


(Al-21) 


(Al-22) 


(Al-23) 


A theorem  due  to  Ricci  states  (which  can  be  proved  quite  simply  by 
using  (Al-21)  - (Al-23))  that  the  covariant  derivatives  of  the  metric  tensor 


are  zero.  That  is 


• $ ‘ 0 


(Al-24) 


Based  on  the  preceding  formulae,  we  state  the  following  results. 


Let  v be  an  arbitrary  vector,  then 


i 1 3 , r i. 

div  v = v . = — (t'g  v ) 

’ Si  35* 


(Al-25) 


Let  T be  an  arbitrary  second  order  tensor.  Then  its  divergence  is 


a vector.  Writing 


b = div  T , 


the  covariant  components  of  b are  given  by 

\ £ (^k>  - ^ 


(Al-26) 
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The  contravariant  components  are  given  by 
b1  = — (/gTlk)  + r ! TrS 

3 C rS 


(Al-27) 


Let  ^ be  a scalar,  then  its  Laplacian  in  general  coordinates  is  given 


• gJk<-4\  - c 

3£J3£  Jk  3C 


(Al-28) 


For  an  aribtrary  vector,  the  following  formulae  establish  the  relations 


between  the  covariant  and  contravariant  components 


Vi  = 8ijV 


i ij 
v = g Jv, 


Expressions  In  Two  Dimensions 


(Al-29) 


In  two  dimensions  both  i and  j vary  from  1 to  2.  We  shall  write 


rl  r r2 

5 = C,  K = n,  = x,  x2  = y 

Denoting  the  partial  derivatives  by  variable  subscripts  we  have 
2 2 

811  = XC  + yC  * 


(Al-30) 


812  = XCXn  + Vn  * 


g22  = xn  + yn 


g = gng22  - (g12)2  = (x4yn  - xnyc)2 
J = , 


. 12  21  , 22  , 

» g ~ 8 “g-io'g  * g = gi  i /g  • 


(Al-31) 


The  Christoffel  symbols  are 

1 ^811  ^811  ^812 

rU  " le22  H + e12<—  "2  ‘ 

r 2 - (g  '-hi  + g (l!22  _2  !hi 

22  t8ll  9n  812 ' 3C  3n 

r 1 = r„  (2  3822x  8g22 

1 22  [g22 3n  3£  ;_g12  3n 


)]/2g 


)]/2g 


2 °12 

rxi  ■ l*u(2-ar 


3gll  3gll 

3n  * -g12  H J/2g 


, 8811  » 3g22. 

^g22  3n  “ g12  3£  ^ 2g 
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2 = r 2 
12  21 

= (gu  - 

1 + r 2 

,ii& 

11 

12 

2g  ac 

■ 1 + r 2 

.-Li* 

12 

22 

2g  an 

81 1 

— — - o — —)/2z 
9£  g12  3n  " 8 


The  Laplacian  of  $ is 

V2*  - lg22*K  - 2g12^n  + «!!♦„„ 

+ ^28i2F12  ” g22ril  " glir22^5 
+ ^2g12ri2  " 822ril  " ®lir22^n^8 


(Al-32) 


(Al-33) 
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Appendix  - A2 

Criteria  For  Elliptic  Equations 


Let  us  consider  two  coupled  nonlinear  partial  differential  equations 


written  generally  with  dependent  variables  x and  y as 
F1(5,n.xIyfx?.xn.y5,ynixK,xen,xn|))  - 0 , 

F2(5,n,x,y,x5,xn,yc,yn,yu,ycn,ynr))  = 0 , 

where  x * x(5,n)  and  y = y(C.n).  For  brevity  we  write 

X4  = Pl’  yC  = p2’  xn  * ql*  yn  ’ q2  * 


(A2-1) 


(A2-2) 


x5£  “ rl*  y££  = r2’  xCn  = Sl’  yCn  = S2  ’ 


x = t,  , y = t0 

nn  l nn  2 

2 2 

If  for  all  pairs  of  real  numbers  (p,v)  such  that  p +v  >0  we  have 


3F1  2 3F1  9F1  2 


(A2-3) 


3F2  2 3F2  3F2 

+3S^PV+3^v2>0  , 


(A2-4) 


then  the  partial  differential  equations  (A2-1)  and  (A2-2)  are  called  elliptic. 
The  Inequalities  (A2-3)  and  (A2-4)  hold  both  for  linear  and  nonlinear  differ- 
ential equations. 

As  an  application  we  consider  Eqs.  (13)  and  (14)  and  apply  (A2-3) 


and  (A2-4) . This  gives 

2 2 
822w  " 2g12uv  + 811V  >0 

Using  the  definitions  of  the  metrical  coefficients  (Al-31),  we  get 

2 2 
(pqj-vpp  + (pq2-vp2)  >0  , 


(A2-5) 


(A2-6) 


which  shows  that  Eqs  (13)  and  (14)  are  uniformly  elliptic  throught  out 


the  domain. 


jtmm i 


Appendix  - A3 

Harmonic,  Subharmonic  and  Superharmonic  Functions 


Let  u(x,y)  be  a continuous  function  in  the  closed  domain  D=DU3D  and 
twice  continuously  differentiable  in  the  open  domain  D.  Consider  the 
equation 


V u = F in  D, 


(A3-1) 


where 


2 2 
v2  „ J_  + Jl 
2 2 
3x 


(A3-2) 


If  F * 0,  then  the  function  u is  called  a harmonic  function  while  if  FiO, 


then  u is  called  a subharmonic  function. 


Maximum  Principle: 


Consider  Eq.  (A3-1)  with  F^O.  From  elementary  calculus  we  know  that 


if  a function  io(x,y)  attains  a maximum  at  (x0»yQ)  then 

to  *0  , a:  *0  , 

x y 

u>  $0  , w ?0  at  (x  ,y  ) , 

xx  * yy  o,Jo'  * 

here  variable  subscripts  denote  partial  differentiation.  Thus,  for  a mixima, 


V2cjf:0 


(A3-3) 


Let  <l>(x,y)  be  a non-negative  function  defined  as 


Defining 


we  find  that 


1 , 2^  2. 

-£t(x  +y  ),  e>0 


U>  = U + <P  , 


V id  *•  F + e }0  , 


since  both  F and  e are  non-negative.  This  contradicts  the  assertion  of  a 
maxima,  i.e.  (A3-3).  Consequently  w can  attain  its  maxima  only  at  the 


boundary  9D.  Also, 


1 2 

(jl)  “ U+$  £ M + -^ep  , 


where  p is  the  radius  of  a circle  containing  D.  Since  u$u>,  hence 


u $ M + ItP 
4 
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Letting  e-*0,  we  find  that  ufM  in  D. 


I 


I 


I 


! 

i 


In  conclusion  we  state  that  if  u satifies  (A3-1)  and  F*0,  then  the 
values  of  u in  D cannot  exceed  the  maximum  on  3D. 

For  harmonic  functions,  viz.,  when  F=o,  we  can  apply  the  above  deduc- 
tions both  to  u and  -u.  This  leads  to  state  that  a harmonic  function  can 
attain  either  its  maximum  or  its  minimum  only  at  the  boundary  3D. 

Subharmonic  Functions: 

The  values  of  a subharmonic  function  u in  a domain  D are  always  below 
the  values  of  the  harmonic  function  v which  coincides  with  u on  3D.  For 
a proof  of  this  assertion  we  follow  Protter  and  Weingerer  [14]. 

Define 

0)  = u - v , 

where  u is  subharmonic,  v is  harmonic,  and  v=u  on  3D.  Thus  to  is  subharmonic 
and  vanishes  on  3D.  According  to  the  maximum  principle  a ><0  in  D.  Thus 

u £ v in  D, 

which  proves  the  assertion. 

Superharmonic  Functions: 

For  the  case  when  F^O  we  similarly  establish  a minimum  principle  which  states 
that  a superharmonic  function  u cannot  attain  a minimum  below  the  values 
of  a harmonic  function  v which  coincides  with  u on  3D. 
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